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ABSTRACT 

We perform a spectrophotometric analysis of galaxies at redshifts z = 4 — 6 in 
cosmological SPH simulations of a A cold dark matter (ACDM) universe. Our models 
include radiative cooling and heating by a uniform UV background, star formation, 
supernova feedback, and a phenomenological model for galactic winds. Analysing a 
series of simulations of varying boxsize and particle number allows us to isolate the 
impact of numerical resolution on our results. Specifically, we determine the luminosity 
functions in B , V, R , i', and z’ filters, and compare the results with observational 
surveys of Lyman break galaxies (LBGs) performed with the Subaru telescope and the 
Hubble Space Telescope. We find that the simulated galaxies have UV colours consistent 
with observations and fall in the expected region of the colour-colour diagrams used 
by the Subaru group. The stellar masses of the most massive galaxies in our largest 
simulation increase their stellar mass from M* ~ 10 11 Mq at * = 6 to M* ~ 10 11,7 Mq 
at z = 3. Assuming a uniform extinction oi E(B — V) = 0.15, we also find reasonable 
agreement between simulations and observations in the space density of UV bright 
galaxies at z = 3 — 6, down to the magnitude limit of each survey. For the same 
moderate extinction level of E(B — V) ~ 0.15, the simulated luminosity functions 
match observational data, but have a steep faint-end slope with a ~ —2.0. We discuss 
the implications of the steep faint-end slope found in the simulations. Our results 
confirm the generic conclusion from earlier numerical studies that UV bright LBGs at 
z > 3 are the most massive galaxies with E(B — V) ~ 0.15 at each epoch. 

Key words: cosmology: theory - galaxies: formation - galaxies: evolution - methods: 
numerical 


1 INTRODUCTION 

Numerical simulations of galaxy formation evolve a comov¬ 
ing volume of the Universe, starting from an initial state 
given by the theory of inflation. Such simulations are in 
principle capable of accurately predicting the properties of 
galaxies that form from these initial conditions, but limited 
computer resources impose severe restrictions on the reso¬ 
lution and volume size that can be reached. In general, it is 
easier to achieve sufficient numerical resolution for high red- 
shift galaxies because the Universe is young and simulations 
are evolved forward in time from the Big Bang. However, 
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observational surveys are often mainly limited to low red- 
shifts, looking outwards and backwards in time from our 
vantage point. In recent years, significant advances in both 
observational and numerical techniques have created an op¬ 
timal overlap range between the two approaches at interme¬ 
diate redshifts (z = 2 — 6), which is therefore a promising 
epoch for comparing theoretical predictions with observa¬ 
tions. This provides a testing ground for the current stan¬ 
dard paradigm of hierarchical galaxy formation in a universe 
dominated by cold dark matter. 

In observational surveys, one of the most important 
techniques for detecting galaxies at redshifts 2 m 3 — 6 makes 
use of the Lyman break, a feature at Any = 4-7Th 3 c/(m e e 4 ) = 
911.7634 A (where m e is the reduced electron mass), the 
wavelength below which the ground state of neutral hydro¬ 
gen may be ionised. Blueward of the Lyman break, a large 
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amount of flux is absorbed by neutral hydrogen, either in the 
galaxy itself or at some redshift along the line of sight. For 
the range we are interested in, the Lyman break is redshifted 
into the optical part of the spectrum. Because of this, these 
galaxies can be detected using optical photometry, making 
them attractive for ground-based surveys; a large difference 
between magnitudes in nearby filters can give an estimate of 
the observer-frame wavelength of the Lyman break, and thus 
the redshift of the galaxy. A galaxy detected in this manner 
is called a Lyman Break Galaxy (LBG). The above ‘break’ 
feature in a redshifted galaxy spectrum causes it to fall in a 
particular location on the colour-colour plane of, e.g., U n — G 
versus G — R colour for 2 sa 3. This colour-selection al¬ 
lows one to preselec t the candidates of high-redshift LBGs 
very efficiently ('e.g.. lSteidel fc Hamiltonll993l : ISteidel et alj 

This method of detecting high-redshift galaxies has 
both advantages and disadvantages. While it is capable of 
detecting a large number of galaxies in a wide field of view 
using relatively little observation time, it cannot assign ex¬ 
act redshifts to galaxies without follow-up spectroscopy. In¬ 
stead, it merely places LBGs into wide redshift bins. More¬ 
over, there is some concern that the procedure may intro¬ 
duce a bias by preferentially selecting galaxies with promi¬ 
nent Lyman breaks. Th ese caveats should be kept in mind 
(e.g., lOuchi et al ] 2QQi when using the results from LBG 
observations for describing the general characteristics of 
galaxies at high redshifts. Nevertheless, the efficiency of se¬ 
lecting high-redshift galaxy candidates coupled with photo¬ 
metric redshift estimates can yield large samples that can¬ 
not be obtained otherwise. Using these techniques as well 
as follow-up spectroscopy, volume li mited surveys of LBGs 
at z > 3 have been constructed (e.g. |LowenthaLet_a ULI99^ 

I Dickinson et a H 120041 : iGiavalisco et alJl200^ L some w ith a 
sample size of N ~ 1000 galaxies fe.g.. ISteidel et al ] 2003 ; 
lOuchi et al.ll200Jl . 

These large datasets make possible interesting compar¬ 
isons with numerical simulations. Perhaps the most impor¬ 
tant fundamental statistical quantity to consider for such a 
comparison is the luminosity function of galaxies; i.e. the dis¬ 
tribution of the number of galaxies with luminosity (or mag¬ 
nitude) per comoving volume. We will focus on this statistic 
here, as well as on the colours, stellar masses, and number 
density of galaxies. 

There have already been several previous stud¬ 
ies of the pr operties of LBGs u sing both semianalytic 
models (e.g. , I Baugh et alJ Il998t ISomerville et alJ l200ll : 
Blaizot_etalJ 120 041) , and cosmologica l simulations (e.g., 
^agamingj200^_[weinber get all l2002l : lHarford fc Gnedira 
2003lkfagaminee^ni2004dl) . As for the semia n alytic mod- 
els of galaxy form ation, both IBaugh et al.1 ill 99Sfl and 
iBlaizot et alJ (I2004T) were able to reproduce the number den¬ 
sity and the correlation function of LBGs, and agree that 
LBGs at 2 ~ 3 ar e massive galaxie s locate d in halos of mass 
~ 10 12 Mq. But ISomerville et al] feoQltl emphasised that 
merger induced starbursts (e.g. Mihos & Hernquist 1996; 
Springel et al. 2005b) could also account for the observed 
number density of LBGs, therefore low mass galaxies with 
stellar masses M* ~ 10 s M@ could also contribute to the 
LBG population. 

Using numerical simulations, iNagamine et a D ll2004dl) 
studied the photometric properties of simulated LBGs at 


z = 3 including luminosity functions, colour-colour and 
colour-magnitude diagrams using the same series of SPH 
simu lations described here (see the erratum of the paper as 
well: INagamine et all2004dl . They found that the simulated 
galaxies have U n — G and G—R colours consistent with obser¬ 
vations (satisfying the colour-selection criteria of Steidel et 
al.), when a moderate dust extinction of E(B — V) = 0.15 is 
assumed locally within the LBGs. In addition, the observed 
properties of LBGs, including their number density, colours 
and luminosity functions, can be explained if LBGs are iden¬ 
tified with the most massive galaxies at z = 3, having typi¬ 
cal stellar masses of M* ~ 10 10 h^o Mq, a conclusion broadly 
consistent with earlier studie s based on hydrodynami c simu¬ 
lations of the ACDM model. INagamine et al 3 •:i2l)05a il: also 
extended the study down to z = 1 — 2 using the same sim¬ 
ulations, and focused on the properties of massive galaxies 

in UV and near-IR wavelengths. 1 _ 

In this paper, we extend the work by INagamine et alJ 
ll2004dh to the LBGs at even higher redshifts z = 4 — 6, fo¬ 
cusing on the colours and luminosity functions of galaxies. 
The paper is organised as follows. In Section |21 we briefly 
describe the simulations used in this paper, and in Sectional 
we outline in detail the methods used to derive the photo¬ 
metric properties of the simulated galaxies. We then present 
colour-colour diagrams in Section 0 discuss stellar masses 
and number densities of galaxies in Section [5] and present 
luminosity functions in Section |S] Finally, we conclude in 
Section |T] 


2 SIMULATIONS 

The simulations analysed in this paper were performed with 
the Smoothed Particle Hydrodynamics (SPH) code Gadget- 
2 (Springel 2005), a Lagrangian approach for modeling hy¬ 
drodynamic flows using particles. W e emp loy the ‘entropy 
formulation’ of SPH JSpringel fc Herntudsl3f2f)0^h which al¬ 
leviates numerical overcooling problems present in other for¬ 
mulations (see e.g. Hernquist 1993; O’Shea et al. 2005). 
The simulation code allows for star formation by converting 
gas into star particles on a characteristic timescale deter¬ 
mined by a subreso lution model for the interstellar medium 
llSpringel fc Hernauistl Il2003 a) , which is invoked for suffi¬ 
ciently dense gas. In this model, the energy from super¬ 
nova explosions adds thermal energy to the hot phase of 
the interstellar medium (ISM) and evaporates cold clouds. 
Galactic winds are introduced as an extension to the model 
and provide a channel for transferring energy and metal- 
enriched mater ial out of the potential wells of galaxies (see 
ISpringel fc Hernauistl I2fl03a . for more details). Finally, a 
uniform U V background radiation field is pr esent, with 
a mo dified lHaardt fe MadatJ \ 19961) spectrum iDave et all 
Il999f ; iKatz et^Lf^96t l. The simulations are based on a 
standard concordance ACDM cosmology with cosmologi¬ 
cal parameters (fi m , Ha, fib, as, h. 70 ) = (0.3, 0.7, 0.04, 0.9, 1), 
where /170 = Ho/(70 km s _1 Mpc -1 ). We assume the same 

1 After the sub mis sion of our manuscript, a preprint by 
iFinlato" et alJ ll2005l) appeared on astro-ph, which studied the 
properties of simulated galaxies at 2 : = 4 in the same simulations 
used in this paper. Overall their results agree well with ours where 
comparison is possible. 
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Size 

-^box 

Run 

N P 

“DM 

fTT'gas 

A i 

N( 6) 

N( 5) 

N( 4) 

N( 3) 

Large 

142.9 

G6 

486 3 

8.99 x 10 s 

1.38 x 10 s 

7.61 

12570 

27812 

48944 

74414 

Medium 

48.2 

D5 

324 3 

1.16 x 10 8 

1.80 x 10 7 

5.96 

6767 

11991 

18897 

24335 

Small 

14.3 

Q6 

486 3 

8.99 x 10 5 

1.38 x 10 5 

1.17 

9806 

12535 

15354 




Q5 

324 3 

3.03 x 10 6 

4.66 x 10 5 

1.76 




7747 


Table 1. Simulation parameters and number of galaxies identified. The simulation boxsize Lt> 0 x gi ven i n Mq . The (initial) number 
of gas particles Np is equal to the number of dark matter particles, so the total particle count is twice Np. The mass of each particle 
i m DM f° r dark matter and m ga s for gas) is given in /iTq 1 A/q . The softening length A£ is given in h^kpc. The last four columns give 
the number of galaxies, N(z), found in the simulation at a given redshift z by the group finder. For the ‘Small’ boxsize, simulation Q6 
was used for z > 4, while Q5 was used for z = 3. 


cosmological parameters when the estimates of the effective 
survey volumes are needed for the observations. 

For this paper, we analyse the outputs of simulations 
with three different box sizes and mass resolutions at red- 
shifts z = 3 — 6. The three simulations employed be- 
l ong to the G-series, D-seri es, and Q-series described in 
ISpringel fc Hernauistl t2 003b ), with corresponding box sizes 
of 142.9, 48.2, and 14.3 h,jQ Mpc in comoving coordinates. 
We shall refer to them as ‘Large’, ‘Medium’, and ‘Small’ 
simulations, respectively. The primary differences between 
the three runs are the size of the simulation box, the number 
of particles in the box, and hence the mass of the individ¬ 
ual particles (i.e. the mass resolution). The parameters for 
each run are summarised in Table 1. The same simulations 
were used for the study of the co s mic star formation his - 
tory llSpringel fc Hernauistl l2003bl: iNa gamine et in 12004 . 
LBGs at z = 3~Eag^ni^^^aiIl2QQ4dr) , d amped Lyman- 
a systems dNagamine. ^8^^^^^fcj^^iauist]|200 4allbll , mas¬ 
sive galaxies at z = 2 l|hhigamin^£^rn]2005a|Ia), and the 
intergalactic medium llFurlanetto et al.ll2004blcidlalf . 

Having three different simulation volumes allows us to 
assess the effect of the boxsize on our results. The measured 
luminosity functions suffer from two different types of res¬ 
olution effects. On the bright end, the boxsize can severely 
limit a proper sampling of rare objects, while on on the faint 
end, the finite mass resolution can prevent faint galaxies 
from being modelled accurately, or such galaxies may even 
be missed entirely. We will discuss these two points more 
more explicitly in Section 5. 


3 METHOD 

Galaxies were extracted fro m the simulation by mean s of a 
group finder as described in iNagamine et alj ( 2004dl) . The 
group finder works by first smoothing the gas and stellar par¬ 
ticles to determine the baryonic density field. The particles 
which pass a certain density threshold for star formation are 
then linked to their nearest neighbour with a higher density, 
unless none of the 32 closest particles has a higher density. 
This is similar to a friends-of-friends algorithm, except it 
does not make use of a fixed linking length. 

After the galaxies are identified, several steps are taken 
to derive their photometric properties. First, we compute the 
spectrum of constituent star particles based on their total 
mass an d metallicity, using a mod ern population synthesis 
model of lBruzual fc Chariot! l2003t) . The spectral energy dis¬ 
tribution (SED) is given as a set of ordered pairs (A, L\(X)). 


Observer-frame wavelength [A] 


0 2000 4000 6000 8000 10000 12000 



Figure 1. Sample LBG spectrum, taken from the ‘Medium’-size 
simulation at z — 4. The bottom axis gives the intrinsic wave¬ 
length of the spectrum, and the top axis is the wavelength red- 
shifted by a factor of 1 + 2 . Also shown are the response functions 
of five Subaru filters (Johnson B, V, and R, and SDSS i' and z'), 
positioned at their effective wavelengths. 


The sampling resolution of this function varies, but in the 
region of interest it is approximately 10A. A sample spec¬ 
trum from the ‘Medium’-size simulation (D5) at 2 = 4 is 
shown in Fig. Q 

Once the intrinsic SED for each source is computed, it 
must be transformed to represent the spectrum as it would 
be seen by an observer on Earth. This process involves sev¬ 
eral steps. Firs t, we app ly the Calzetti dust extinction law to 
the spectrum fcalzetti et alJEoOO i. which accounts for in¬ 
trinsic extinction within the galaxy. The specific values we 
adopt for the strength of the extinction will be discussed in 
more detail below. Next, w e redsh ift the spectrum and ac¬ 
count for IGM absorption (Madaui (1995). Finally, we com¬ 
pute the photometric magnitudes by convolving the result¬ 
ing SED with various filter functions. This allows us to deter¬ 
mine the apparent magnitude of each object for commonly 
employed filters in the real observations. 

The formula used for this computation may be derived 
as follows. For a given source (defined by its flux per unit 
frequency /„(^)), observed through a given filter (defined by 
its filter r esponse function RM), the apparent magnitude is 
given by dFukugita et al.lfl995l . Eq. 7): 


m = —2.5 log 


f f v (v)R(v)d\nv 

J Cv(v)R{v) dlni/’ 


(1) 


where C v {v) is the reference SED. For the AB magnitude 
system, C v is a constant (10 -19 ' 44 erg s -1 cm -2 Hz -1 ), and 
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for the Vega system, C v is the SED of the star Vega. The 
above formula may be rewritten in terms of wavelength using 
the relation f v (u) = (A 2 /c)/a(A), where c is the speed of 
light, and the observed flux may be related to the intrinsic 
luminosity by /a (A) = La(A /(1 + z))/[47rd|(l + 2)], where 
A is the wavelength in the observer frame, and <1l is the 
luminosity distance to redshift 2 . This gives (for cgs units): 


m.AB = —2.5 log 


/ A^a R(X) dA 
47rd| c(l + 2 ) f jR( A) dA 


48.60, 


( 2 ) 


The absolute magnitude Mab may be determined from this 
equation by setting 2 = 0 and d,L = 10 pc. For monochro¬ 
matic magnitudes, as shown in Fig.Q the equation reduces 
to Mab{ A) = — 2.5log(A 2 Z/A) + 13.83, assuming A is given 
in units of [A] and L\ is in units of [Lq A 1 ]. 

Absorption by dust and extinction by the IGM each 
add a multiplicative factor to f\ as a function of wave¬ 
length inside the integral. For dust absorption, the factor is 
lQ-fc(A/(i+ 2 ))B^ w j lere A:(A/(1 + 2 )) is the Calzetti extinction 
function, and E = E(B—V) is the extinction in B—V colour, 
taken to be a free parameter. There is no simple theoretical 
constraint on E(B — V) except that it must be nonnega¬ 
tive, so we simply consider a range of values for E(B — V) 
to study the ex tinction effect system atically. Since the lat¬ 
est surveys ('e.g. lShanlev et al.ll200ih suggest that E(B — V) 
ranges from 0.0 to 0.3 with a mean of ~ 0.15, we adopt three 
fiducial values of E(B — V ) = 0.0, 0.15, and 0.30. In most of 
our figures, they will be indicated by the colours blue, green, 
and red, respectively. We discuss a different choice for the 
assignment of extinction to galaxies in Section [^] 

For the IGM extinction, the factor is exp[—r(A, 2 )], 
where t(A, 2 ) is the effec tive optical depth owing to both 
contin uum dMadaulll995 l footnote 3) and line extinction 
dMadaulll995l Eo. 15): 

r{X,z) = 0.25xi' 46 (a 0A6 — 1 ) 

+ a:*' 68 (9.4a 018 + 0.7a -132 - 0.023a 1 ' 68 - 10.077) 

+ E A ^^) 346 ’ ( 3 ) 

3 =2 J 


where x c = X/Xl v , a = (1 + z)/x c , and the Aj are the line 
strength coefficients. We consider only the four strongest 
lines, corresponding to j = 2 to 5. This absorption becomes 
highly significant in the blue bands at redshifts greater than 
3, as shown in Fig. 0 

Note that dust absorption is applied in the rest frame of 
the galaxy, while IGM extinction is applied in the observer’s 
frame. The integration can also be done in the rest frame 
rather than the observer frame by substituting A with A(1 + 
2 ). Thus the overall formula to compute m from L\ is: 

mAB = -48.60 - 2.5 log ( 1 

\ 4™ £ c 

/ ALa(A) 10- fc(A)B e -r(>( 1 +*).*) R( A(1 + 2 )) dA\ 

/ i R(X(1 + 2 )) dA ) ^ 


We have used several filters for our calculations, each 
defined by a response function R( A). Specifically, for com¬ 
parison with observations from the Subaru telescope, we 
used their filters B, V, R, i !, and 2 ' (Johnson -M organ - 
Cousins system; see Section 2.6 of iMivazaki et a DSqqI), 
which provide good coverage of all optical wavelengths, and 



Figure 2. Total optical depth owing to IGM absorption, at red- 
shifts of 6, 5, 4, and 3, taking into account continuum absorption 
and absorption from the four strongest lines. The right axis gives 
the corresponding increase in apparent magnitude. The names 
of the five Subaru filters are positioned at their effective wave¬ 
lengths. 


some into the near-infrared. lOuchi et aL d2()11 1 i treated the 
i' magnitude as the standard UV magnitude for 2 ^ 4, and 
the z' magnitude as the standard UV magnitude for 2 ~ 5. 

For comparison with surveys that did not use the Sub¬ 
aru filters, and for a more general UV luminosity function, 
we used a boxcar-shaped filter (i.e. response function set 
equal to unity) centered at 1700 A and with a half-width of 
300A, in the rest frame of the observed galaxy. Note that this 
is actually a different filter in the observer’s frame depend¬ 
ing on the redshift of the observed object (e.g., for a 2 = 4 
object, it is centered at rest frame 8500A, and for a 2 = 5 
object, it is centered at 10200A). We refer to the magnitude 
measured with this filter simply as the ‘UV-magnitude’ in 
this paper. Depending on a particular survey’s capabilities, 
observationally determined UV magnitudes may be based on 
a slightly different wavelength than 1700A, but it is a rea¬ 
sonable assumption that the resulting magnitudes are com¬ 
parable. 


4 COLOUR-COLOUR DIAGRAMS 

LBGs are identified based on a significantly dimmer mag¬ 
nitude in a filter blueward of their Lyman break compared 
with a filter redward of their Lyman break. This difference 
is manifest as a significantly redder colour. Moreover, two 
filters redward of the Lyman break should not exhibit ab¬ 
normal dropouts with respect to one another, a fact that 
can distinguish them from interlopers with very red spec¬ 
tra. Thus, in order to select LBGs from the sample, colour 
selection criteria are very important. For instance, galaxies 
at 2 ~ 4 will have Lyman breaks at approximately 4600A, 
between the B and R filters (as shown in Fig. Q. so these 
galaxies will have large B — R colours. But, we also expect 
them to have moderate R—i' colours, since both R and i! are 
redward of 4600 A. The exact colour-colour selection criteria 
used by each survey are determined empirically by placing 
a sample of spectroscopically identified LBGs on a colour- 
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Figure 3. Colour-colour diagram for the plane of B — R versus R — i' . The selection criteria for B.Rz-LBGs in the Subaru survey are 
shown by the dashed lines. The contours in panels (a) and (b) show simulated galaxies in the ‘Large’ (G6) run for z — 3 and z = 4, 
respectively. Blue, green, and red contours correspond to E(B — V) = 0.0, 0.15, and 0.30, respectively. 
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Figure 4. Colour-colour diagram for the plane of V — i' versus i' — z'. The selection criteria for Viz- l.BGs in the Subaru survey are 
shown by the dashed lines. The contours in panels (a) and (b) show the simulated galaxies in the ‘Large’ (G6) run for z = 4 and 2 = 5, 
respectively. Blue, green, and red contours correspond to E(B — V) = 0.0, 0.15, and 0.30, respectively. 


colour diagram, or computing the track of local galaxies with 
known spectra (or artificial spectra of galaxies generated by 
a population synthesis model with assumed star formation 
histories) as a func tion of redshift. 

For example, lOuchi et alJ ( 2004 *1 use the following 
colour criteria for selecting 2 ~ 4 galaxies: 

B-R > 1.2 (5) 

R-i! < 0.7 (6) 


B-R > 1.6(7? — i') + 1.9 (7) 

Galaxies identified as LBGs by these criteria are called BRi- 
LBGs. Similar criteria exist for V — *’ versus i' — z' colours, 
and R — i! versus i! — z' colours; galaxies selected in this way 
are called Viz- LBGs and Tfe-LBGs, respectively. Each of 
these three classes of LBGs corresponds to an approximate 
range in redshift space. The central redshifts for BRi, Viz, 
and Riz populations are approximately 4, 5, and 5, respec- 
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Figure 5. Colour-colour diagram for the plane of R — i' versus i' — z'. The selection criteria for Riz LBGs in the Subaru survey are 
shown by the dashed lines. The contours in panel (a) and (b) show the simulated galaxies in the ‘Large’ (G6) run for z = 4 and 2 = 5, 
respectively. Blue, green, and red contours correspond to E(B — V) = 0.0, 0.15, and 0.30, respectively. 



miB (i'-band) 

Figure 6. Colour-magnitude diagram for the ‘Medium’ boxsize 
(D5) run at z = 4. Blue, green, and red contours outline the 
distribution of galaxies for E(B — V) = 0.0, 0.15, and 0.30, re¬ 
spectively. Also shown in green are vertical lines approximately 
delineating the number of star particles in galaxies. Most galax¬ 
ies found to the left of the first line comprise 1600 star particles 
or more, and 200 particles mark the end of our interval of confi¬ 
dence. Black dashed lines are cutoffs for BRi LBGs as observed 
by the Subaru group. The vertical line indicates the 3 a limiting 
magnitude of i' = 26.9, and the horizontal line represents one of 
the colour-colour selection criteria used to identify LBGs. 


tively. The Riz selection has a narrower range of redshift 
than the Viz selection, so we will use the BRi sample to 
compare to our 2 = 4 simulations, and Riz sample to com¬ 
pare to our 2 = 5 simulations. 

In Figs. E] |21 we show the colour-colour diagrams of 
our simulated galaxies. Overall, the agreement with the ob¬ 
servation is good, and the simulated galaxies fall within the 


same regi on as the observed gala xies, a conclusion consis¬ 
tent with lNagamine et alJ 1 21III Id :. Very few 2 = 3 galaxies 
in our simulation would be detected as BRi- LBGs, but a 
large fraction of 2 = 4 galaxies would. Similarly, relatively 
few 2 = 4 galaxies in our simulations would be detected as 
V* 2 -LBGs or J?i 2 -LBGs compared with 2 = 5 galaxies. This 
result appears to be relatively insensitive to the amount of 
Calzetti extinction, at least for the range of extinction values 
we considered. 

In Fig. 0 we also show the colour-magnitude diagram 
on the plane of i'-band apparent magnitude and R—i! colour 
for the ‘Medium’ (D5) run. This figure shows that all the 
simulated galaxies brighter than uiab (*' — band) = 27 sat¬ 
isfy the colour-selection of R — i! < 0.7. Since the brightest 
galaxies in the simulations are the most massive ones, this 
means that the LBGs in the simulations are the brightest 
and most massive galaxies with E(B — V) ~ 0.15 at each 
epoch. The situation of course changes when a larger value 
of extinction is allowed, as such galaxies could become red¬ 
der than R — i' = 0.7. Such dusty starburst galaxies may 
exist in the real universe, but we are not considering them 
in this paper by restricting ourselves to E(B — V) < 0.3. 


5 GALAXY STELLAR MASSES 

Fig ,|7|shows stellar mass M* vs. UV magnitude for the simu¬ 
lated galaxies over several redshifts and for different extinc¬ 
tion values, in relation to the survey limiting magnitudes. 
From the figure it is seen that the ‘Large’ box size simula¬ 
tion contains galaxies as massive as Af* ~ 10 11 /i)T 0 1 Mq at 
2 = 6, and A/* ~ lO 11,1 at 2 = 3. Larger objects are 

too rare to be found in a simulation of this size. The diago¬ 
nal lines in the figure depicting mass-to-light ratio show that 
this value is generally increasing going from higher to lower 
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Figure 7. Contour plots of stellar mass of galaxies vs. UV magnitude, for redshifts z = 3 — 6, using the ‘Large’ (solid colour contour) box 
size. Blue, green, and r ed contours represen t extinctions of E(B — V) = 0.0, 0.15, and 0.3 respectively. The dashed black lines indicate 
the magnitude limits of lSteidel et alJ J2003il at z = 3, the SDF sample for z — 4 and 5, and the HST GOODS for z = 6. Diagonal lines 
show lines of constant stellar mass to light ratio; the value of / (\L\) is labeled for each line, in units of Mq/Lq. 


redshift, so that the luminosity per stellar mass is decreasing 
with time. 

Fig. El shows the cumulative number density (integral 
of the luminosity function for galaxies brighter than a cer¬ 
tain magnitude limit) of galaxies at various redshifts and for 
different extinction values for the ‘Large’ (solid curves) and 
‘Medium’ (dashed curves) box sizes. Every simulation of lim¬ 
ited size will underpredict the continuum value by a certain 
amount, depending on the bright-end cutoff imposed by the 
finite volume, so it is not surprising that the ‘Medium’ sim¬ 
ulation gives systematically lower densities than the ‘Large’ 
simulation. For c o mpari son, we sho w values dete r mined by 
lAdelberger et all l2Q03|) at 2 = 3- jOuchi et alJ (l2004l at 
z = 4 and 5, and lBoiiwens et, alJ (l2004fl at z = 6. The best 


match with the observations is reached for E(B — V) = 0.15 
at 2 = 3, but the required extinction appears to slightly 
increase towards E(B — V) = 0.3 at higher redshifts. 

At z = 3 and for E(B — V) = 0.15, we find a value for 
the number density of n(R < 25.5) ~ 1 x 10~ 3 (/iRMpc)~ 3 . 
This magnitude was determined by the limiting magnitude 
in the survey of lAdelberger et alJ f 2003 1: for values at dif¬ 
ferent magnitudes, refer to Fig. El Similarly, at z = 4 and 5, 
we determine values of n(i' < 26.5) ~ 6 x 10 _3 (/iRMpc) -3 
and n(z' < 26.0) ~ 1.5 x 10 _3 (/iRMpc) _3 . The variation in 
these resul ts prima rily reflects the different limiting magni¬ 
tudes of lOuchi et alJ (' 2004 1 instead of an internal evolution 
of the LF over this redshift range. Finally, at z = 6, we mea¬ 
sure a value of n(z' < 29.0) ~ 2.1 x 10~ 2 (hRMpc) -3 , where 
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Figure 8. Total number density of galaxies below a threshold magnitude in the ‘Large’ box size (G6 run, solid curves) and ‘Medium’ 
box size (D5 run, dashed curves) simulations. Blue, green, and red curves represent the extinction values E( B — V) = 0.0, 0.15, and 0.3, 
respectively. The vertical solid black lines roughly represent the magnitude limits of the lSteidel et all J2003i) sample for z = 3, the SDF 
sample for z = 4 and 5, and the HST GOODS for z = 6. 


the li miting magnitude was chosen as in iBouwens et alJ 
<2004 . 


6 LUMINOSITY FUNCTIONS 

The most widely used analytic parameterisation of the 
galaxy lu minos ity function is the Schechter Function 
fechechteill 1974 . the logarithm of which is given by 

log($(M)) = log(0.41n(10)$*) + fj,(a + 1) - 10"/ ln(10), (8) 

where /x = — 0.4(M — M*), and $*, M*, and a are the nor¬ 
malisation, characteristic magnitude, and faint-end slope, 
respectively. We note that throughout this section, we plot 
luminosity functions (LFs) in terms of magnitude rather 
than luminosity. Brighter objects will thus appear farther 
left on the abscissa than fainter objects. 


The LF measured from our simulations suffers both 
from boxsize and resolution effects, so we expect it to be 
physically meaningful only for a certain limited range of 
luminosities. At the faint end, objects are made up by a 
relatively small number of particles, and may not be well- 
resolved, and even smaller objects will be lost entirely. In 
our simulations, the luminosity functions generally have a 
peak at around 100 stellar particles. The turn-over on the 
dim side of this peak owes to the mass resolution. In order to 
avoid being strongly affected by this limitation, we usually 
discard results based on galaxies with fewer than 200 stellar 
particles. 

At the bright end, objects become increasingly rare. 
When there is only of order one object per bin in the entire 
box, the statistical error of the LF dominates and we can¬ 
not reliably estimate the abundance. In order to improve the 
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Figure 9. Typical luminosity function, depicting our approxi¬ 
mate criterion for physical relevance of the data, and uncertainty. 
The top and right axes show simulation-based measurements of 
objects: the approximate number of stellar particles which make 
up an object in the given magnitude bin, and the number of 
objects in the entire box in a given magnitude bin (of size 0.2 
mags). The bottom and left axes are the corresponding physical 
measurements: the magnitude as derived from the spectrum of 
the objects, and the number density as a luminosity function. A 
fiducial cutoff of 200 stellar particles is chosen as the lower limit 
for accurately resolved objects, although this value is somewhat 
arbitrary. The dotted line indicates Poisson V~N uncertainties, 
and the dashed line indicates the standard deviation of the mean 
between the luminosity function derived for eight disjoint sub- 
regions of the entire box. Since the latter uncertainty is greater 
throughout most of the range, only it is shown on subsequent 
plots. Both forms of uncertainty are dependent on the bin size, 
0.2 mags. 


sampling of these objects, a larger simulation volume needs 
to be chosen, which is however in conflict (for a given par¬ 
ticle number) with the usual desire to obtain a good mass 
resolution. 

Fig. m illustrates these numerical limitations and de¬ 
fines a region where the LF results can be trusted. The cut¬ 
offs we indicated here are however approximate. A precise 
determination of the interval over which the measured LF 
is physically significant requires running many simulations 
with successively higher resolution, and looking for an in- 
terva l of convergence. S uch a p rogramme was carried out 
bv iSpringel fc Hernauist (2003bj) in their study of the cos¬ 
mic star formation rate. We here analyse three different box 
sizes taken from their set of simulations, yielding a good 
coverage of magnitude space. Where appropriate, we also 
combine these results to obtain a measurement covering a 
larger dynamic range. In subsequent figures, the approxi¬ 
mate interval of physical significance for the LF measure¬ 
ments will be shown with a thicker line than the rest of the 
curve. Strictly speaking, it is only this interval that can be 
compared reliably with observations. 

Also shown in Fig. Elis the estimated uncertainty in the 
LF owing to sample size, which we estimated in two ways. 



E(B-V) 


Figure 10. Correlation between UV magnitude and extinction 
value E(B — V). Blue points, plo tted with respect t o the left axis, 
show the original data from |2haEleY_et_alJ 12001 *1. Green lines 
show the best linear fit to the data, and the lines denote 1 -a scat¬ 
ter. Black contours outline the extinction values we assigned to 
the simulated galaxies using the procedure described in Section IH1 
(for the ‘Large’ boxsize G6-run at z = 4) in order to match up 
with this correlation. The red histogram, plotted with respect to 
the right axis, shows the distribution of extinction values of the 
black contours. However, for the reasons described in the text, we 
do not use this method of assigning variable E(B — V) hereafter. 


First, we calculated vV Poisson statistical errors, simply 
by taking the square root of the number of objects in each 
bin. Second, we divided the box into eight octants, and com¬ 
puted the standard deviation of the mean between the LFs 
computed with each of the eight octants. In all cases, the sec¬ 
ond method produced larger uncertainty over the interval of 
physical significance, indicating that the ‘cosmic variance’ 
error owing to our limited boxsize exceeds a simple Poisson 
estimate. We therefore use error estimates obtained with the 
octant method in our subsequent figures on the LF results 
and ignore the Poisson errors. 

Typically, we assumed three values for the extinction, 
E(B — V ) = 0.0, 0.15, and 0.30, and produced LF-estimates 
for them separately. In this procedure, we hence always as¬ 
signed a single extinction value to every galaxy for which 
we computed magnitudes. However, in the real universe 
we instead expect a distribution of e xtinction values (e.g., 
IShaolev et all200ltlQuchi et all2004i . which could be quite 
broad. This prompted us to explore possible effects owing to 
a ‘variable extinction’. To this end, we first introduced ran¬ 
dom scatter into the values for extinction: Instead of appling 
the same extinction to all galaxies, each galaxy was assigned 
an individual value of E(B — V) determined by a Gaussian 
random variable with a mean of 0.15 and a standard de¬ 
viation of 0.10. Moreover, a cutoff was imposed so that no 
galaxy had a negative extinction value. We found that such 
variable scatter tends to smooth out the luminosity func¬ 
tion somewhat, as expected, but it does not produce results 
readily distinguishable from a uniform extinction value. 

There is evidence for a correlation between UV mag - 
nitude and extinction. In particular, IShaolev et a.l] ( 2001 1 
found that, over the magnitude range studied, dust obscures 
almost exactly enough flux to give all galaxies a similar 
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Figure 11. i'-band luminosity functio n for both ‘Large’ (G6, panel a) and ‘Medium’ boxsize (D5, panel b) simulation runs at z = 4, and 
the observational data for BRi -LBGs iOuchi et al .12004 Fig. 16). Simulation LFs are plotted as blue, green, and red curves representing 
extinction values of E(B — V) = 0.0, 0.15, an d 0.30, respectively. The Subaru survey data are shown as black crosses and boxes for the 
two different survey fields iOuchi et all2004 Fig. 16). 



Figure 12. Luminosity function for all three boxsizes at z = 4. 
Extinction colour-coding is the same as in the previous plots. 
Solid and dashed lines are best Schechter fits to Subaru BRi- 
LBG data, assuming a value for th e faint -end slope of a = —2.2 
and -1.6, respectively JOuchT et al • 20041 Fig. 16). The dashed- 
dotted line has a slope of -2.0, but is not a fit to the Subaru 
data. 

apparent magnitude. The data from IShanlev et~aT. ( 1200 1). 
along with the extinction values we adopted to emulate it, 
appear in Fig. ITUl When these extinction values are used, 
brighter galaxies have larger values of E(B — V). The ef¬ 
fect of this distribution of extinction values is such that the 
brightest galaxies become dim enough to agree with the ob¬ 
served magnitudes, while the fainter ones suffer little extinc¬ 
tion so that they agree with the narrow range of observed 


magnitude. Taken together, this effect produces a simulated 
LF that matches the empirical one quite well within the ob¬ 
served range of magnitudes, requiring only moderate extinc¬ 
tion at the faint end. However, invoking variable extinction 
in this manner is of course bound to succeed at some level, 
because we here ‘hide’ the difference between simulated and 
observed LFs in the variable extinction law. While such a 
law in principle may exist, we prefer here to systematically 
study the effect of extinction by assuming different values of 
E(B — V) uniformly for the entire sample. 

In Figs. mi through m we show the LFs derived from 
the simulated observations as coloured curves with data 
points from observational surveys overplotted. These are the 
main results of this paper. We first compare the i'-band 
LFs directly to the observations. For z = 4, this is shown in 
Fig. llll for the ‘Large’ (G6) and ‘Medium’ boxsize (D5) runs, 
along with data from lOuchi et all f 20041) . The observational 
data points do not extend faint enough to offer any overlap 
with the ‘Small’ (Q5 & Q6) run’s LF, therefore we do not 
show the results from the ‘Small’ simulation in this figure. 
As seen in the figures, the data points are roughly consis¬ 
tent with the green curves, corresponding to an extinction 
value of E(B — V) ~ 0.15. The result of our ‘Large’ run cov¬ 
ers the entire range of the observed magnitude ranges, but 
the ‘Medium’ run overlaps only with the fainter side of the 
Subaru data because of the moderate simulation volume. 

We then compared the simulated LFs for all three box- 
sizes (‘Large’, ‘Medium’, and ‘Small’) to the best-fitting 
Schechter function as determined from the Subaru survey. 
For z = 4, this is shown in Fig. The two Schechter func¬ 
tions plotted in this figure correspon d to the best fit to the 
Subaru data obtained by lOuchi et al.1 (120041) with a fixed 
faint-end slope of either a = —1.6 (dashed line) or —2.2 
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Figure 13. z '-band luminosity function for both ‘Large’ (G6, panel a) and ‘Medium’ boxsize (D5, panel b) simulations at z = 5. Survey 
results for similar redshifts are shown with symbols. Simulation data are plotted as blue, green, and red curves representing extinction 
values of E(B — V) = 0.0, 0.15, and 0.30, re spectively. Ouchi et al.’s js'-band survey data are shown as black crosses and boxes for the 
two different survey fields iOuchi et all2004 Fig. 16). The updated I -band survey data of llwata et al 2 Hoi are shown as triangles. 



Figure 14. ^'-band luminosity function for all three boxsizes at 
z = 5. Extinction colour-coding is the same as in the previous 
plots. Solid and dashed lines are best fit Schechter functions to 
Subaru Riz -LBG data, assuming a value for t he faint-end slope 
of a = —2.2 and —1.6, respectively iOuchi et al • 200j Fig. 16). 
The dashed-dotted line has a slope of —2.0, but is not a fit to the 
data. 

(solid line). The simulated LF with E(B — V) = 0.15 is 
consistent with the faint-end slope between these two val¬ 
ues, and the value of a = —2.0 (dash-dotted line) appears 
reasonable. 

Similar results emerged for z = 5. In Fig. we show 
the ‘Large’ an d ‘Medium’ boxsize LFs and compare them 
with data from IOuchi et, al] ll2004l and llwata et alJ ll2004l . 
At 2 = 5, the survey data appear to be more consistent with 


the simulations using E(B — V) = 0.3, particularly at the 
bright end, suggesting a slightly larger extinction at 2 = 5 
than at z = 4. And again, when all three boxsizes are plot¬ 
ted along with the best-fitting Schechter function in Fig. rm 
the most consistent value for the faint end slope a appears 
to lie between —1.6 and —2.2. Note that there are sm all dis- 
crep ancies between the LF-estimates bv lOuchi et alJ (120041 
and llwata et alJ (l2004l . This may owe to slight differences 
in the colour-selection criteria used in the two surveys. 

In Fig.H3 we compare the UV-magnitude LFs at z = 6 
in th e ‘Large’ an d ‘Medium’ boxsize simulations to data from 
iBouwens et alJ 1120041 . Again, a value of E(B — V) between 
0.15 and 0.30 leads to the best match with the data, and an 
extinction with E(B — V) = 0.15 is favoured based on this 
comparison. 

Finally, Fig. 1161 examines the evolution of the LF over 
the redshift range studied. To this end, we plot all redshifts 
for all boxsizes using the single extinction value E(B — V) = 
0.15. Overall, the LF of our simulations shows little if any 
evolution over the redshift range in question. The absence 
of strong evolution is probably related to the fact that the 
evolution of the cosmic star formation rate (SFR) density 
is quite mil d from z = 3 to 6 i n our s imulations, as dis¬ 
cussed by ISpringel fc Hernauistl 12003b ) and Hernquist & 
Springel (2003). In both SPH and total variati on diminish¬ 
ing (TVD) simulations llNagamine et a hi 12004 1 the cosmic 
SFR continu es to rise gradually from 2 = 3 to 5, and peaks 
at 2 = 5 — 6. iNagamine et al.l f 2005afl have shown that the 
evolution of the LF from 2 = 3 to 2 = 2 is about 0.5 mag, so 
it is perhaps not too surprising that the evolution at higher 
redshift is of comparably small size. Note that in terms of 
proper time, the redshift interval from 2 = 6 to 3 is only 
about as long as the interval from 2 = 3 to 2. 
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Figure 15 . UV lum inosity function at z = 6 for both ‘Large’ and ‘Medium’ boxsize simulations and the observational data from 
iBouwens et alJj2004l black crosses). Simulation data are plotted as blue, green, and red curves, representing extinction values E(B—V) = 
0.0, 0.15, and 0.30, respectively. 



Figure 16. Rest frame UV luminosity function for all three box- 
sizes for redshifts z = 6 (red), 5 (orange), 4 (green), and 3 (blue). 
For clarity, results for each LF are only plotted over the reliable 
range, and error bars are omitted. The extinction value for all 
curves is here E(B — V) = 0.15. 

7 DISCUSSION & CONCLUSIONS 

Using state-of-the-art cosmological SPH simulations, we 
have derived the colours and luminosity functions of sim¬ 
ulated high-redshift galaxies and compared them with ob¬ 
servations. In particular, we have employed a series of simu¬ 
lations with different boxsizes and resolution to identify the 
effects of numerical limitations. We find that the colours 
of galaxies at z = 4—6 agree with the observed ones 
on the colour-colour planes used in observational studies, 
and our results con firm the generic conclusion from earlier 
numerical studies dNagaminel 120021 : IWeinberg et alJ l2002t 


iNagamine et al . 2004dl . 2005al ) that UV bright LBGs at 
z > 3 are the most massive galaxies with E(B — V) ~ 0.15 
at each epoch. 

The simulated LFs are in good agreement with the data 
provided an extinction of E(B—V) = 0.15 — 0.30 is assumed. 
The faint-end slope of our results is consistent with a value 
between a = —1.6 and —2.2, as found with the Subaru data 
llOuchi et alJ 120041) . The simulated LFs best match a very 
steep faint-end slope of a ~ —2.0. The recent analysis of 
Hubble Space Telescope Ultra Deep Field (UDF) data by 
lYan fc Windhorst 1 2004 ) also suggests a quite steep faint- 
end slope of a = —1.8 to —1.9, in good agreement with the 
simulations. 

The steep faint-end slopes at z > 6 found here are in¬ 
teresting because they have significant implications for the 
history of the reionisation of t he Un iv erse. As discussed 
by, e. g., ISonnge^^Hem^uis 3 ( 20031 1). lYan fc Windhorstl 
(2004:), and lCen et al.l (12005^ 1. it may be possible to reionise 
the Universe at 2 = 6 with ionising photons from Population 
II stars in normal galaxies alone if the faint-end slope of the 
galaxy lu minosity function is sufficiently steep («< — 1.6). 
Note that iNagamine et a D d2004dl) also found a similarly 
steep faint-end slope at z = 3 in the same SPH simulations. 

The suggestion that the faint-end slope may be much 
steeper at high redshifts than in the Local Universe is in¬ 
triguing, and the agreement between the recent observations 
and our simulations lends encouraging support for this pro¬ 
posal. However, this immediately invokes the question of 
what mechanisms changed the LF over time and gave it 
the shallow slope (a ~ —1-2) observ ed locally in the 2dF 
llCole et al.ll200ll) and SDSS surveys llBlanton et, ahlEoPlI) . 
A simil arly strong flattening d oes not occur in our simu¬ 
lations llNagamine et al.ll2004dl) . but this could owe to an 
incomplete modeling of the physics of feedback processes 
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from supernova and quasars. We note that the SPH simula¬ 
tions employed in this paper already inclu ded a galactic wind 
model fsee ISnringel fc Hernauistll2003al . for details) which 
drives some gas out of low-mass halos, but the effect is not 
sufficiently sensitive to galaxy size to produce a significantly 
flattened faint-end slope. The shallowness of the faint-end 
slope of the local LF hence remains a challenge for cosmo¬ 
logical simulations and may point to the need for o ther p hys¬ 
ical processes, such as black hole growth ('e.g-.ISpringel e t alJ 
l2005bl : iDi Matteo et, al ]1 200 -ot iRobertson et, alfl^OORlh 

In light of this, it is therefore particularly encourag¬ 
ing that the simulations are more successful at high red- 
shift. The steep faint-end slope found by observations in 
this regime is an important constraint on the nature of the 
feedback processes themselves. Presently, the observational 
results still bear substantial uncertainties, however, stem¬ 
ming mostly from their limited survey volume which makes 
them prone to cosmic variance errors, and less from the faint¬ 
ness levels reached, which already probe down to itlab = 29 
magnitude in the case of the UDF survey. The next gener¬ 
ation of wide and deep surveys will shed more light on this 
question in the near future. 
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